﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace 高斯正反算
{
    public class Guass
    {
        public int m_zone_index;
        private int m_zone_wide;
        private double m_lat_origin;
        private double m_e0;
        private double m_e1;
        private double m_e2;
        private double m_e3;
        private double m_e;
        private double m_es;
        private double m_esp;
        private double m_ml0;
        protected double m_central_meridian;
        protected double m_false_east = 500000;
        protected double m_false_north;
        protected double m_scale_factor;
        protected double m_semi_major = 6378137.0;
        protected double m_semi_minor = 6356752.3142451793;

        private readonly double radToDegree = 0.017453292519943295;
        private readonly double degreeToRad = 57.295779513082323;

        public Guass()
        {
            this.m_zone_wide = 3;
            this.m_zone_index = -1;
            this.m_scale_factor = 1.0;
            this.m_lat_origin = 0.0;
        }
        public bool setCentralLon(int zone_width, int zone_index)
        {
            if (zone_width != 3)
            {
                if (zone_width != 6)
                {
                    return false;
                }
                this.m_central_meridian = (double)(zone_width * zone_index - 3) * radToDegree;
            }
            else
            {
                this.m_central_meridian = (double)(zone_width * zone_index) * radToDegree;
            }
            this.m_zone_index = zone_index;
            this.m_zone_wide = zone_width;
            return true;
        }
        public bool setZoneWide(int wide)
        {
            this.m_zone_wide = wide;
            return true;
        }
        public bool GeoToProject(double longitude, double latitude, ref double coord_x, ref double coord_y)
        {
            coord_x = 0.0;
            coord_y = 0.0;
            double num = radToDegree * longitude;
            double num2 = radToDegree * latitude;
            this.m_es = 1.0 - Math.Pow(this.m_semi_minor / this.m_semi_major, 2.0);
            this.m_e = Math.Sqrt(this.m_es);
            this.m_e0 = 1.0 - 0.25 * this.m_es * (1.0 + this.m_es / 16.0 * (3.0 + 1.25 * this.m_es));
            this.m_e1 = 0.375 * this.m_es * (1.0 + 0.25 * this.m_es * (1.0 + 0.46875 * this.m_es));
            this.m_e2 = 0.05859375 * this.m_es * this.m_es * (1.0 + 0.75 * this.m_es);
            this.m_e3 = this.m_es * this.m_es * this.m_es * 0.011393229166666666;
            this.m_ml0 = this.m_semi_major * (this.m_e0 * this.m_lat_origin - this.m_e1 * Math.Sin(2.0 * this.m_lat_origin) + this.m_e2 * Math.Sin(4.0 * this.m_lat_origin) - this.m_e3 * Math.Sin(6.0 * this.m_lat_origin));
            this.m_esp = this.m_es / (1.0 - this.m_es);
            if (this.m_zone_wide == 6)
            {
                int num3 = (int)(longitude / (double)this.m_zone_wide) + 1;
                this.m_central_meridian = radToDegree * (double)(num3 * this.m_zone_wide - (this.m_zone_wide >> 1));
            }
            else
            {
                int num3 = (int)((longitude - 1.5) / (double)this.m_zone_wide) + 1;
                this.m_central_meridian = radToDegree * (double)num3 * (double)this.m_zone_wide;
            }
            double num4 = this.adjust_lon(num - this.m_central_meridian);
            double x = Math.Sin(num2);
            double num5 = Math.Cos(num2);
            double num6 = num5 * num4;
            double num7 = Math.Pow(num6, 2.0);
            double num8 = this.m_esp * Math.Pow(num5, 2.0);
            double num9 = Math.Tan(num2);
            double num10 = Math.Pow(num9, 2.0);
            double d = 1.0 - this.m_es * Math.Pow(x, 2.0);
            double num11 = this.m_semi_major / Math.Sqrt(d);
            double num12 = this.m_semi_major * (this.m_e0 * num2 - this.m_e1 * Math.Sin(2.0 * num2) + this.m_e2 * Math.Sin(4.0 * num2) - this.m_e3 * Math.Sin(6.0 * num2));
            coord_x = this.m_scale_factor * num11 * num6 * (1.0 + num7 / 6.0 * (1.0 - num10 + num8 + num7 / 20.0 * (5.0 - 18.0 * num10 + Math.Pow(num10, 2.0) + 72.0 * num8 - 58.0 * this.m_esp))) + this.m_false_east;
            coord_y = this.m_scale_factor * (num12 - this.m_ml0 + num11 * num9 * (num7 * (0.5 + num7 / 24.0 * (5.0 - num10 + 9.0 * num8 + 4.0 * Math.Pow(num8, 2.0) + num7 / 30.0 * (61.0 - 58.0 * num10 + Math.Pow(num10, 2.0) + 600.0 * num8 - 330.0 * this.m_esp))))) + this.m_false_north;
            return true;
        }
        public bool ProjectToGeo(double coord_x, double coord_y, ref double longitude, ref double latitude)
        {
            longitude = 0.0;
            latitude = 0.0;
            this.m_es = 1.0 - Math.Pow(this.m_semi_minor / this.m_semi_major, 2.0);
            this.m_e = Math.Sqrt(this.m_es);
            this.m_e0 = 1.0 - 0.25 * this.m_es * (1.0 + this.m_es / 16.0 * (3.0 + 1.25 * this.m_es));
            this.m_e1 = 0.375 * this.m_es * (1.0 + 0.25 * this.m_es * (1.0 + 0.46875 * this.m_es));
            this.m_e2 = 0.05859375 * this.m_es * this.m_es * (1.0 + 0.75 * this.m_es);
            this.m_e3 = this.m_es * this.m_es * this.m_es * 0.011393229166666666;
            this.m_ml0 = this.m_semi_major * (this.m_e0 * this.m_lat_origin - this.m_e1 * Math.Sin(2.0 * this.m_lat_origin) + this.m_e2 * Math.Sin(4.0 * this.m_lat_origin) - this.m_e3 * Math.Sin(6.0 * this.m_lat_origin));
            this.m_esp = this.m_es / (1.0 - this.m_es);
            int num = (int)(coord_x / 1000000.0);
            if (this.m_zone_wide == 6)
            {
                this.m_central_meridian = radToDegree * (double)(num * this.m_zone_wide - this.m_zone_wide / 2);
            }
            else
            {
                this.m_central_meridian = radToDegree * (double)num * (double)this.m_zone_wide;
            }
            double num2 = coord_x - this.m_false_east;
            double num3 = coord_y - this.m_false_north;
            double num4 = (this.m_ml0 + num3 / this.m_scale_factor) / this.m_semi_major;
            double num5 = num4;
            long num6 = 0L;
            while (true)
            {
                double num7 = (num4 + this.m_e1 * Math.Sin(2.0 * num5) - this.m_e2 * Math.Sin(4.0 * num5) + this.m_e3 * Math.Sin(6.0 * num5)) / this.m_e0 - num5;
                num5 += num7;
                if (Math.Abs(num7) <= 1E-10)
                {
                    break;
                }
                num6 += 1L;
            }
            if (Math.Abs(num5) < 1.5707963267948966)
            {
                double x = Math.Sin(num5);
                double num8 = Math.Cos(num5);
                double num9 = Math.Tan(num5);
                double num10 = this.m_esp * Math.Pow(num8, 2.0);
                double num11 = Math.Pow(num10, 2.0);
                double num12 = Math.Pow(num9, 2.0);
                double num13 = Math.Pow(num12, 2.0);
                num4 = 1.0 - this.m_es * Math.Pow(x, 2.0);
                double num14 = this.m_semi_major / Math.Sqrt(num4);
                double num15 = num14 * (1.0 - this.m_es) / num4;
                double num16 = num2 / (num14 * this.m_scale_factor);
                double num17 = Math.Pow(num16, 2.0);
                double num18 = num5 - num14 * num9 * num17 / num15 * (0.5 - num17 / 24.0 * (5.0 + 3.0 * num12 + 10.0 * num10 - 4.0 * num11 - 9.0 * this.m_esp - num17 / 30.0 * (61.0 + 90.0 * num12 + 298.0 * num10 + 45.0 * num13 - 252.0 * this.m_esp - 3.0 * num11)));
                double num19 = this.adjust_lon(this.m_central_meridian + num16 * (1.0 - num17 / 6.0 * (1.0 + 2.0 * num12 + num10 - num17 / 20.0 * (5.0 - 2.0 * num10 + 28.0 * num12 - 3.0 * num11 + 8.0 * this.m_esp + 24.0 * num13))) / num8);
                longitude = degreeToRad * num19;
                latitude = degreeToRad * num18;
            }
            else
            {
                longitude = degreeToRad * (3.1415926535897931 * (double)((num3 < 0.0) ? -1 : 1) / 2.0);
                latitude = degreeToRad * this.m_central_meridian;
            }
            return true;
        }
        private double adjust_lon(double x)
        {
            long num = 0L;
            double num2 = 2147483647.0;
            double num3 = 4.61168601E+18;
            while (Math.Abs(x) > 3.1415926535897931)
            {
                if ((int)Math.Abs(x / 3.1415926535897931) < 2)
                {
                    x -= (double)((x < 0.0) ? -1 : 1) * 3.1415926535897931 * 2.0;
                }
                else
                {
                    if ((double)((int)Math.Abs(x / 3.1415926535897931 * 2.0)) < num2)
                    {
                        x -= (double)((int)(x / 3.1415926535897931 * 2.0)) * 3.1415926535897931 * 2.0;
                    }
                    else
                    {
                        if ((double)((int)Math.Abs(x / (num2 * 3.1415926535897931 * 2.0))) < num2)
                        {
                            x -= (double)((int)(x / (num2 * 3.1415926535897931 * 2.0))) * (6.2831853071795862 * num2);
                        }
                        else
                        {
                            if ((double)((int)Math.Abs(x / (num3 * 3.1415926535897931 * 2.0))) < num2)
                            {
                                x -= (double)((int)(x / (num3 * 3.1415926535897931 * 2.0))) * (6.2831853071795862 * num3);
                            }
                            else
                            {
                                x -= (double)((x < 0.0) ? -1 : 1) * 3.1415926535897931 * 2.0;
                            }
                        }
                    }
                }
                num += 1L;
                if (num > 4L)
                {
                    break;
                }
            }
            return x;
        }
    }
}
